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In this study we analyze the phase and group velocity of three-dimensional linear traveling waves 
in two sheared flows, the plane channel and the wake flows. This was carried out by varying the 
wave number over a large interval of values at a given Reynolds number inside the ranges 20 — 100, 
1000 — 8000, for the wake and channel flow, respectively. Evidence is given about the possible 
presence of both dispersive and non-dispersive effects which are associated with the long and short 
ranges of wavelength.We solved the Orr-Sommerfeld and Squire eigenvalue problem and observed 
the least stable mode. It is evident that, at low wave numbers, the least stable eigenmodes in the 
left branch of the spectrum beave in a dispersive manner. By contrast, if the wavenumber is above 
a specific threshold, a sharp dispersive to non-dispersive transition can be observed. Beyond this 
transition, the dominant mode belongs to the right branch of the spectrum. The transient behavior 
of the phase velocity of small three-dimensional traveling waves was also considered. Having chosen 
the initial conditions, we then show that the shape of the transient highly depends on the transition 
wavelength threshold value. We show that the phase velocty can oscillate with a frequency which 
is equal to the frequency width of the eigenvalue spectrum. Furthermore, evidence of intermediate 
self-similarity is given for the perturbation field. 


I. INTRODUCTION 

The relationship between the oscillation frequency and 
the wave vector - which is known as the dispersion rela¬ 
tion - is a general relationship that is important in all 
areas of wave physics. Simply stated, a dispersion rela¬ 
tion is the function a;(fc) for a harmonic wave. For the 
most basic waves, where the speed of propagation c is 
a constant, the dispersion relation is simply uj{k) = ck. 
That is, the angular frequency is a linear function of the 
wave vector. The ratio u/k is the propagation speed c 
and it is known as the phase velocity. In this case, the 
phase speed is independent of k. However, this is not 
always the case. An example is given by the propagation 
of light in a dielectric medium, where the index of refrac¬ 
tion n = Co jc (where cq is the speed of light in a vacuum) 
depends upon the wavelength. Another case of an inter¬ 
esting, nonlinear dispersion relation is found in the wave 
function that describes a free particle propagating in a 
given direction with a given momentum. 

In the context of temporal linear perturbation dynam¬ 
ics of viscous incompressible flows, the general solution to 
the Orr-Sommerfeld equation at a fixed Reynolds number 
yields a dispersion relation between a real wavenumber k 
and a complex angular frequency to. Any arbitrary initial 
disturbance can be decomposed and expressed as a linear 
combination of eigenfunctions. For bounded flows, the 
Orr-Sommerfeld spectrum consists of an infinite number 
of discrete eigenvalues and it is complete [1]. While for 
unbounded flows, only a finite number of discrete eigen¬ 
values exists and thus a continuum of modes must exist 
[2, 3]. 


* Accepted for publication on Physical Review E 
Email address for correspondence: daniela.tordella@polito.it 


In the long term the least stable mode (LSM), i.e. the 
mode which has the maximum temporal growth rate, pre¬ 
vails. In temporal stability analysis, the focus is usually 
placed on the range of wavenumbers where the maximum 
temporal growth rate is positive and in general the in¬ 
formation on the relation uii = u}i{k) is more abundant 
with respect to the corresponding ujr = uij.(k) [2, 4-8]. 

Spectral theory of hydrodynamic waves has been de¬ 
veloped to a high degree of sophistication, with intri¬ 
cate analytical tools supplemented by accurate multi¬ 
dimensional computational methods [9-13]. Neverthe¬ 
less, unresolved key issues in this classical problem can 
be still identified. For instance, let us consider the fol¬ 
lowing questions: 

- Perturbations at a fixed wavelength and Reynolds num¬ 
ber don’t always retain the same frequency of oscillation 
during their transient. Thus the phase speed may vary 
inside transients. How are these variations linked to the 
structure of the eigenvalue spectra, in particular to their 
branched structure? 

- In open sheared flows, like wakes for instance, at fixed 
values of the control parameter, differences in the disper¬ 
siveness associated to the set of eigenvalues belonging to 
the left and right branches of the spectrum are observed. 
A fact which has not been yet examined in the case of 
wall flows. In synthesis, what happens to the least stable 
mode (LSM) when the perturbation wavelength is varied 
by fine steps over a large range? Does the dispersive¬ 
ness change when the LSM moves from one branch to 
the other? 

- Is intermediate self-similarity of the perturbation pro¬ 
files present inside transients? 

Specihcally, we have considered these issues to take aim 
at two points. Firstly, the relationship between the an¬ 
gular frequency corresponding to the least stable Orr- 
Sommerfeld mode and the perturbation wavenumber was 
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analyzed, in a wide range of wavenumbers, that is from 
k = 0.20 to fc = 20. This interval extends from below 
to well above the eventual instability range. We con¬ 
sider two classical incompressible shear flows, the plane 
channel and wake flows. The analysis was carried out by 
considering Reynolds numbers in the range 1000 — 8000 
for the channel flows, and in the range 20 — 100 for the 
wake. The aim was to obtain both the phase speed and 
the group velocity and by comparison observe possible 
variations of the dispersive behavior with the wavenum¬ 
ber. 

We observed the existence of a wavenumber thresh¬ 
old {kd, in the following) that separated the waves 
which propagate in a dispersive way from the waves 
that propagate in a non-dispersive way. While of inter¬ 
est for the study of wave packets, the existence of this 
dispersive/non-dispersive transition can also affect the 
transient dynamics of the traveling waves. 

Secondly, the temporal evolution of sheared flow small 
perturbations is investigated, with a particular empha¬ 
sis on the phase velocity transients and self-similarity is¬ 
sues. Indeed, the frequency or phase velocity transient 
has been poorly investigated to date. For instance, in the 
wake flow, attention has mainly focused on the frequency 
of vortex shedding for the most unstable spatial scales 
[14-16]. Only very recently, subcritical wake regimes 
(up to values 30% below the critical value) of transiently 
amplified perturbations have been studied by consider¬ 
ing the spatio-temporal evolution of wave packets [17]. 
The situation is quite different within the context of at¬ 
mosphere and climate dynamics. Here, the interaction 
between low-frequency and high-frequency phenomena, 
which is related to the existence of very different spatial 
and temporal scales, is believed to be one of the main 
reasons for planetary-scale instabilities [18, 19]. However, 
due to the inherent strong nonlinearity, the evolution of 
single scales cannot be observed in geophysical systems 
and thus also these studies usually do not account for the 
structure of the frequency transient of a single wave. In 
this work, the temporal evolution between the early tran¬ 
sient and the temporal asymptotic state was considered 
in detail. It is observed that this transition can be char¬ 
acterized by phase velocity jumps inside the perturbation 
temporal evolution. This implies that the perturbation 
may experience acceleration or deceleration. This be¬ 
havior depends on the initial condition adopted and on 
the wavelength of the perturbation, in particular if this 
is smaller or greater with respect to the dispersive/non- 
dispersive threshold. 

The organization of the paper is as follows. In section 2 
both an eigenvalue problem and an initial-value problem 
are formulated, (see also Appendix A). The dispersion 
relation of the least stable mode is discussed in Section 
3. The phase velocity behavior in the transient is then 
described in Section 4. In Section 5, a discussion on 
the existence of the intermediate term and its empirical 
determination is presented. Conclusion remarks follow in 
Section 6. 



FIG. 1. Basic flows and perturbation scheme in the x-y phys¬ 
ical plane. A Cartesian reference frame is adopted, with unit 
vectors ei, e 2 , es in the x, y, z directions, respectively. The 
base flow profiles are qualitatively represented in black with 
arrows. A generic perturbation wave is represented in red. 
k = oei -I- yes is the wavenumber vector, 0 is its angle with 
respect to the basic flow \J = U(y)ei. In the physical domain, 
the longitudinal space variable is unbounded for the channel 
flow and positive for the wake flow. However, due to the 
adoption of the near-parallel flow approximation, see section 
II, the domain for the plane wake is also seen as unbounded. 
The spanwise coordinate is unbounded for both flows. Note 
that the transversal computational domain width is defined 
through the parameter y/. For the channel flow yj = ±1, 
while for the wake flow i// is defined so that the numerical 
solution be insensitive to further extensions of the computa¬ 
tional domain {yf = 20 for short waves and j// up to 100 for 
long waves). 


II. PHYSICAL PROBLEM AND 
FORMULATION 


We consider two typical shear flows, the plane chan¬ 
nel flow, an archetype of bounded flows, and the plane 
bluff-body wake, one of the few free flow archetypes. A 
Cartesian coordinate system is adopted, with origin at 
the channel mid plane in the first case and at the bluff- 
body location, for the wake case. The x, y, z axis are ori¬ 
ented in the streamwise, the transversal and the spanwise 
directions, respectively (see Fig. 1). After introducing 
arbitrary small perturbations the linearized, viscous and 
incompressible governing equations in non-dimensional 
form read: 
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dxU + dyV + dzW = 0 ( 1 ) 

dfU + U dxU + vU'+ dxP = ( 2 ) 

Re 

dtv + U dxv + dyp = ( 3 ) 

He 

dtw+ UdxW + dzP = ( 4 ) 


where (u{x,y,z,t), v{x,y,z,t), w{x,y,z,t)) and 

p{x,y,z,t) are the components of the perturbation 
velocity and pressure, respectively. The domain is 
—oo < x,z < oo, — 1 < y < 1 for the channel flow, and 
0 < a: < 00 , —oo < z < oo, —oo < ?/ < oo for the wake. 
No-slip boundary conditions for the channel flow imply 
vanishing perturbation at the walls, u = v = w = Qss, 
2 / = ± 1 , while decaying disturbances are considered for 
the wake flow, which is unbounded, so that u,v,w ^ 0 
as 2 / —>■ oo. Periodicity boundary conditions are used 
in both X and z directions. The channel half-width h, 
and the body diameter D are considered as reference 
external length scales. The reference velocity for the 
channel flow is the centreline velocity C/cl, while in 
the wake case the free-stream velocity C// is considered. 
The reference time is the convective one. Consequently, 
the flow control parameter is the Reynolds number, 
defined as Re = Uci^h/u, for the channel flow, and as 
Re = UfD/iy for the wake flow, where v is the kinematic 
viscosity. The channel base flow is represented by the 
plane Poiseuille solution 

U{y) = l-y\ (5) 

As a wake basic flow, we use the first two order terms 
of the Navier-Stokes solution described in [20, 21] and 
reported below: 

U{y- a:o. Re) = Co - (g) 

where Cq = 1 and Ci = 1.22 -h O.OOOOGTRe^. The 
two-dimensional wake is a spatially evolving free flow. 
Leaving aside the near field, that is highly non paral¬ 
lel since it hosts the two symmetric countercirculating 
vortices that constitute the separation region, the inter¬ 
mediate and long term wake is a near-parallel flow. The 
wake slowly becomes thicker according to a law which, at 
first order, scales as {x/Re)^'^. As representation of this 
steady subcritical flow we consider the previously men¬ 
tioned asymptotic expansion solution in inverse powers 
of X which was also used in convective instability wake 
intermediate asymptotics analysys [22-24]. In particular, 
we consider the intermediate - far field well represented 
by sections located in the interval x G [5, oo], in the range 
of Reynolds numbers [20,100]. The basic flow is approxi¬ 
mated by freezing it at three longitudinal stations placed 
at xq = 10, 20, 50. In so doing, the basic flow is parame¬ 
terized with the downstream station xq and the Reynolds 


number. It is thus homogeneous in x and z. The mo¬ 
mentum equations Eq. (2-4) can be expressed in terms 
of the velocity-vorticity formulation: 

[(dt + udx)v^-u"dx-^v^]v = o (7) 

[dt + Udx - = -U'dzV (8) 

uJy = dzU-dxW ( 9 ) 

Eqs. (7,8,9) are solved by means of a combined Fourier- 
Fourier (channel) and Laplace-Fourier (wake) transform 
in the plane normal to the basic flow profile. The trans¬ 
formation of a generic variable 'g{x, y, z, t) reads: 

+ 00 +00 

5(2/,/;a,7) = y J g{x,y,z,t)e~^°‘^~'‘^''dxdz (10) 

— OO 0, —oo 

where a and 7 are the streamwise and spanwise 
wavenumbers, respectively. For the sake of simplicity, in 
the present work only real wavenumbers are considered, 
even if for the wake flow the formulation allows to take 
into account spatially evolving waves. The wavenumber 
modulus is fc = -|- 7 ^ and the wave obliquity angle 

is (f) = tan~^(y/a), see Fig. 1. The governing equations 
in the wavenumber space are thus formulated: 

[(5t -b iaU) {dl - k^) - iaU” 

= ^ ( 11 ) 

[(a, + ^aU) - {dl - k^)]iOy = -i^U'v (12) 

The streamwise and the spanwise velocity components, 
^ can be recov¬ 

ered from the continuity equation and from the vorticity 
definition. The boundary conditions associated to the 
system 11-12 for the channel flow are 

i)(±l, t) = dyv{±.l, t) = a)y (±1, t) = 0, (13) 

while in the wake case we consider hnite-energy harmonic 
velocity as \y\ -G 00 , and vanishing vorticity in the free- 
stream, see also [2] and [25]: 

dlv = k'^v, Ldy=0 y ^ ±oo,yt (14) 

About the initial conditions, please, refer to section IV. 
The eigenvalues a = ar + ici and the eigenfunctions of 
the eigenproblem associated to the system Eq. (11-12) 
are computed by means of three different methods: a 
5*^-order Galerkin method based on Chandrasekhar ex¬ 
pansions (see [26], and Appendix A) , a finite difference 
4*^-order scheme, [27], and a hybrid spectral collocation 
method based on Chebyshev polynomials [ 8 ] . A compar¬ 
ison of this techniques and the details of the first method 
are given in Appendix A and B. 
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Two different numerical methods are used to solve 
the initial value problem: the method of lines, based on 
finite difference spatial discretization and a Runge-Kutta 
(2,3) temporal integration scheme, and the 5*^-order 
Galerkin method based on the Chandrasekhar functions 
expansion cited above. The results obtained by the 
two approaches present a really good agreement. The 
computational domain is [—?//,j//]. For the channel 
flow yf = 1 , while for the wake flow y/ is defined so 
that the numerical solution be insensitive to further 
extensions of the computational domain size {yf = 20 
for short waves and y/ up to 100 for long waves). The 
solution computed by the initial-value formulation at 
long enough time shows excellent agreement with the 
behavior predicted by the modal analysis. 


In the case of the wake we use as reference transversal 
observation point yo = 1 or yo = 5, and in the case of 
the channel flow the point yo = 0.5. The frequency [25] 
is thus 

w(t;yo,Q;,7) = \dd{t;yo,a,j)\/dt. (20) 

The phase velocity is defined as 

c = (a;/fc)k, (21) 

where k = cos{(j))ei + sin((())e 3 is the unit vector in the 
direction of k. 


In order to measure the growth of the perturbations, 
we define the kinetic energy density as 

1 f+Vf 

e{f,a,j) = — {\uf + \if + \^\^)dy. (15) 

4y/ J-vf 

We then introduce the amplification factor, G, as the ki¬ 
netic energy density normalized with respect to its initial 
value. 


G(t;a, 7 ) = e{t;a,j)/e{t = 0 ;a, 7 ). (16) 

Since the temporal asymptotic behavior of the linear per¬ 
turbations is exponential, the temporal growth rate r is 
defined [7] as 

r{t; a, 7 ) = log{G)/{2t). (17) 

The temporal evolution of the kinetic energy is given 
by the equation [ 8 , 12 ]: 

de 
dt 

+ \dyi^y\+k^\l^y\'^^dy, ( 18 ) 

where the bar indicates the complex-conjugate and 3 is 
the imaginary part. While the dissipative term is always 
negative, the convective term has undetermined sign and 
can be responsible for kinetic energy growth. 

The frequency of the perturbation uj is defined as 
the temporal derivative of the unwrapped wave phase 
0 (y,t;a, 7 ), at a specific spatial point along the y direc¬ 
tion. The wrapped phase. 


J U' (^vdyV - jvujy^ dy 


dw{y,t-,a,j) = arg{v{y,t-,a,-f)), (19) 

is a discontinuous function of t defined in [—tt, -I-tt], while 
the unwrapped phase 6 * is a continuous function obtained 
by introducing a sequence of 27r shifts on the phase val¬ 
ues in correspondence to the periodical discontinuities. 


III. DISPERSIVE TO NON-DISPERSIVE 
BEHAVIOR IN THE LONG TERM 

The focus of this section is on the long-term temporal 
behavior of small amplitude traveling waves, and on the 
relation between the frequency of the least stable eigen¬ 
value and the wavenumber. The wavenumber is within a 
wide range 0.2 < k < 20, with a small step Sk = 0.005. 
This allows to highlight a transition from dispersive to 
non-dispersive behavior. This transition is also useful to 
understand the temporal evolution of general arbitrary 
initial perturbations, which is the object of the following 
sections. The results here presented, see Tables I and II 
and Fig. 2, have been compared with literature data ob¬ 
tained by means of different methods of investigation, in 
particular laboratory, modal and Initial Value Problem 
analysis. Least stable eigenvalue frequencies are deduced 
from the computation of the eigenvalue spectra at con¬ 
stant wavenumber, waveangle and Reynolds number. 
Figure 2 (a,b) shows the dependence of the phase veloc¬ 
ity c on the polar wavenumber k in the range k € [0.2, 5]. 
In both the flows we observe the existence of a thresh¬ 
old wavenumber that we indicate with kd where a step 
variation of the phase velocity occurs. For k > kd the 
phase velocity is approximately constant and is nearly 
equal to the group velocity, a fact which highlights a 
non-dispersive behavior. For k < kd instead, c depends 
on A: in a different way from the group velocity and the 
behavior is dispersive. 

In other word, inside the range of wavelengths that can 
be hosted in the system, there is a wavenumber thresh¬ 
old above which the least stable mode is always on the 
tip of the right branch of the eigenvalue spectrum (about 
the branched structure of the eigenvalue spectrum, see 
for example [7, 8 ] and panels d and h of figure 4 be¬ 
low). By increasing the wavenumber, the frequency of 
the least stable mode increases in a proportional way, 
which yields non-dispersiveness. Viceversa, below this 
threshold, the least stable mode is located on the tip of 
the left branch. By varying the wavenumber one sees that 
this least stable mode frequency does not vary in a way 
that is proportional to the wavenumber variation, which 
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FIG. 2. Dispersion relations for longitudinal waves {(j> = 0) in both the channel and wake flow. Results from eigenvalue 
analysis, in particular from least stable mode characteristics. (a,b) Phase and group velocity as a function of the wavenumber 
and comparison with other studies. The phase velocity is represented with the blue line, while the group velocity line is red. 
For more info on the relation tu(fc) see also Figs. SI and S2 in the SM. The wavenumber range [0.2,5] is discretized with 
Ak = 0.005. kd indicates the threshold value that separates the dispersive waves from non-dispersive waves, (a) Channel flow 
with Re = 5000, kd = 1.62 and comparison with other studies, see Grosh & Salwen 1968 [4], Ito 1974 [28], Nihsioka et al. 1975 
[29], Asai & Florian 2006 [30]. (b) Wake flow. Re = 50 and xo = 20. Here kd = 1.32. A comparison is given with the phase 
velocity values found by Roshko 1974 ]31], Nihsioka & Sato 1974 [32], Norberg 1994 [33], Williamson 1989 [14], Paranthoen 
et al. 1999 [15], Pier 2002 [34], Barkley 2006 [35], Belan & Tordella 2006 [22], Tordella et al. 2006 [23], Giannetti & Luchini 
2007 [36]. Wherever not available, the wavelength has been set equal to the mean of the values measured by Paranthoen et 
al. (A = 7D) and Williamson (A = 7.15D). (c-f) Convective and dissipative energy variation terms as dehned in eq. 22 as 
a function of the wavenumber, for the two least-stable eigenmodes of the left and right branch of the spectrum. Blue lines 
represent the least-stable solution. For the channel flow at Re = 5000, the evolution of two dispersive, one of which centered 
on kd, and one non-dispersive wave packets is shown in three movies in the SM online [37]. 


yields dispersivity. In terms of phase and group velocity 
this threshold yields a sharp transition because the tips 
of the left and right branches of the eigenvalue spectrum 
have substantially different frequencies, see again Fig. 4, 
panels d,h. It should be considered that during the tran¬ 
sient part of the evolution of the sum of a number of 
perturbations with different wave numbers (for example 
a wave packet) one can better perceive the contribution 
of the least stable modes for each wavenumber, which in 


turn can be associated to different spectral branches and 
have different frequencies and growth or decay factors. In 
this last case, in the long term, the solution will mainly 
contain the less stable components for each specific wave 
number inside. These least stable modes will disperse if 
their wavelength is below the kd threshold slightly sensi¬ 
tive to the Reynolds number, see tables I and II, and that 
the opposite is true for wavelengths above the threshold. 
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A question naturally arises: why does the set of least- 
damped waves become non-dispersive above a given fc? 
In this regard, we like to highlight the following. As 
already said, it can be observed that for k < kd, the 
least-stable Orr-Sommerfeld eigenvalue is located on the 
tip of the left branch of the eigenvalues spectrum. As 
described in further detail in Section V (see figures 6 and 
7), the related velocity eigenfunction is a shear mode: 
the highest perturbation variations (dv/dy) are located 
where the base flow vorticity dU/dy is high (in the shear 
region). This region is located close to the walls for the 
plane Poiseuille flow, while it is close to the inflection 
points in the wake flow. We refer to these perturbations 
as wall modes for channel flow, and in-wake modes for the 
wake. Vicevera, for k > k^, the least-stable eigenvalue 
is located on the tip of the right branch and the domi¬ 
nant perturbation solution is an external mode: in this 
case the perturbation varies rapidly where base flow has 
low vorticity (out of the high-shear region). In synthesis, 
when the wavelength of a perturbation is short enough 
{k > kd) the dominant mode variation dv/dy is confined 
in a region of low base flow vorticity. In this case, note 
that the dispersion relation for the least-damped modes 
is close to the dispersion relation for a uniform base flow 
which is always non-dispersive (see [8], chap 3.2). Here, 
we refer to central modes for the channel flow, and out- 
of-wake modes for the wake flow. 

Coming back to the comparison with data in litera¬ 
ture, one can see that the agreement between labora¬ 
tory data as well as numerical modal analysis and our 
work is very good. In particular, as regards the chan¬ 
nel flow, we observe a good agreement with the works 
of Grosh & Salwen 1968 [4], Ito 1974 [28], Nishioka et 
al. 1975 [29], and Asai & Floryan 2006 [30] (see panel 
(a) of Fig. 2). Grosh & Salwen give the least stable 
eigenvalues in the range of a [1 - 2.4] , (^ = 0, and in 
the range [5 - 25000] of the Reynolds number (under our 
normalization). Asai and Floryan investigated experi¬ 
mentally the effects of wall corrugation on the stability 
of wall-bounded shear flows and found a good agreement 
with the Orr-Sommerfeld theory coupled to roughness in 
the form of a single Fourier harmonic and in the form of 
spanwise grooves with rectangular and triangular shapes 
[38]. 

For the wake flow, it should be observed that in litera¬ 
ture data are mostly focused on the value of the shedding 
frequency. Our references include the works by Roshko 
1954 [31], Nishioka & Sato 1974 [32], Williamson 1989 
[14], Norberg 1994 [33], Paranthoen et al. 1999 [15], Pier 
2002 [34], Barkley 2006 [35], Belan & Tordella 2006 [22], 
Tordella et al. 2006 [23], and Giannetti & Luchini 2007 
[36]. The agreement is good also in this case. Indeed 
the relative error with respect to the measurements by 
Williamson and Paranthoen et al. is about 3.3%. We 
would like to remind the reader that wake data in litera¬ 
ture usually give information on the frequency at which 
vortex shedding takes place but don’t give information 
about the wavenumber of the shedding. Thus, to ad- 


TABLE I. Values of the dispersive regime threshold 
wavenumber kd for the channel flow, for different Reynolds 
number and obliquity angles. The uncertainty on kd due to 
the discretization is ±0.005. 


Re 

(/) = 0 

Channel flow 
(p ~ 7r/6 

</) = 7r/4 

cf} = 7r/3 

1000 

2.071 

2.111 

2.168 

2.256 

2000 

1.883 

1.922 

1.979 

2.073 

3000 

1.764 

1.803 

1.866 

1.960 

4000 

1.686 

1.725 

1.784 

1.878 

5000 

1.623 

1.662 

1.721 

1.815 

6000 

1.576 

1.615 

1.670 

1.765 

7000 

1.536 

1.568 

1.627 

1.720 

8000 

1.497 

1.536 

1.589 

1.682 


dress this lack of information, we decided to associate to 
the frequency values reported in literature the average 
of the wavenumber values measurements currently avail¬ 
able for Re = 50, that is Williamson [14] and Paranthoen 
[15]. So that, leaving aside the two points related to re¬ 
sults from near-parallel multiple scale convective analysis 
[22, 23], where wavelength and frequency are determined 
along the streamwise direction, all the points in figure 2 
b have been superposed on the same wavenumber. 

Goming back to Table 1 and 11, on can see that the 
threshold kd is a function of the Reynolds number and 
the wave angle but for the wake flow - which is a weakly 
streamwise evolving flow - it is also a function of the 
streamwise station xq- In particular, for the channel 
flow, values of kd have been computed in the range 
Re G [1000,8000] and (f G [0,7r/3]. We observe that 
kd decreases as Re and increases as the wave angle. The 
kd — Re trend is opposite in the wake case, see Table 11 
where for Re G [20,100] the dependence on the flow sta¬ 
tion xq is shown {xq = [10, 20, 50]). It is observed that in 
the case of the channel flow, the system tends to become 
non-dispersive at all wavenumbers as Re tends to infinite. 
As regards to the wake, this trend is reversed. This may 
not be intuitive but one can consider that with increasing 
Re, even though the wake narrows, the shear intensity in¬ 
creases and the range of dispersive waves widens. This 
supports the fact that the base vorticity intensity plays 
a key role in the dispersion. 

Finally, by considering the kinetic energy equation for 
two-dimensional perturbations, 

ft = (“'’»»«) ‘‘y 

" -V-" 

Convection 

~R^ j + + (22) 

"-V-" 

Dissipation 

we compared the convective and the dissipative terms for 
the least stable mode on the left and right branches of 
the spectrum respectively, see Fig. 2 panels (c,d,e,f). For 
simplicity, we have considered only symmetric modes but 
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TABLE II. Values of the dispersive regime threshold 
wavenumber kd (uncertainty of ±0.005) for the wake flow, for 
different Reynolds number, obliquity angles and streamwise 
station xq. 


Re 

o 

II 

-e- 

Xq = 10 

(f) — 7r/6 

= 7r/4 

= 7r/3 

20 

1.073 

1.032 

0.977 

0.868 

30 

1.373 

1.333 

1.276 

1.154 

40 

1.636 

1.598 

1.534 

1.407 

50 

1.875 

1.835 

1.772 

1.640 

60 

2.091 

2.052 

1.989 

1.862 

70 

2.293 

2.253 

2.195 

2.068 

80 

2.485 

2.448 

2.390 

2.264 

90 

2.663 

2.628 

2.575 

2.453 

100 

2.837 

2.802 

2.749 

2.633 

Re 

(/) = 0 

*0 = 20 

(f) — 7r/6 

= 7r/4 

<f) = tt/S 

20 

0.758 

0.735 

0.696 

0.617 

30 

0.974 

0.946 

0.903 

0.817 

40 

1.158 

1.131 

1.087 

0.997 

50 

1.326 

1.298 

1.251 

1.161 

60 

1.478 

1.451 

1.408 

1.318 

70 

1.623 

1.596 

1.553 

1.463 

80 

1.756 

1.733 

1.689 

1.604 

90 

1.885 

1.858 

1.819 

1.737 

100 

2.001 

1.983 

1.944 

1.862 

Re 

(/) = 0 

Xq = 50 

cf) — 7r/6 

= 7r/4 

= 7r/3 

20 

0.482 

0.461 

0.441 

0.401 

30 

0.615 

0.593 

0.573 

0.522 

40 

0.732 

0.714 

0.684 

0.633 

50 

0.840 

0.815 

0.795 

0.734 

60 

0.935 

0.916 

0.886 

0.835 

70 

1.026 

1.007 

0.977 

0.926 

80 

1.112 

1.088 

1.068 

1.007 

90 

1.189 

1.169 

1.148 

1.098 

100 

1.267 

1.249 

1.229 

1.179 


it can be shown that the results do not change when an¬ 
tisymmetric modes are considered. In Fig. 2, Cl and 
D1 indicate, for the least stable mode, the contribution 
of the convective and viscous terms of equation (22) to 
the normalized kinetic energy temporal rate of change, 
respectively. C2 and D2 indicate the contribution given 
by the second-last least damped mode, which for k > kd 
belongs to the P-family of eigenvalues for the Poiseuille 
flow, and to the continuous branch for the wake flow. 
It can be observed that by means of an interchange of 
modes taking place at k = kd, the system, beyond this 
threshold, settles to the mode that minimizes the nor¬ 
malized kinetic energy loss. In synthesis, beyond kd the 
traveling waves are stable and non-dispersive, and the 
rate of kinetic energy loss associated to the convective 
effects is independent (wake case) or very mildly depen¬ 
dent (channel case) on the wavenumber. Cases of mode 
exchange for shear or stratified flows with complex dis¬ 
persion relationship have been reported and analyzed by 
other authors (see, for instance, [3, 39] for the Poiseuille 


pipe flow and and [40] for the stratified Couette flow). 
For a review, see [6], sections 2.2 and 2.6. It should 
be noted that the mode exchange observed here for the 
plane channel flow went mostly unnoticed up to now. In 
fact, we were able to only find a short reference inside 
Fig 3.8b in [8], page 71. Possibly, because it takes place 
in the very stable region above the neutral curve of the 
fc-Re plane, that is, a region not of great interest when 
the focus is on instability. 

In the next sections the transient evolution of small 
perturbations is considered. We show that the threshold 
kd has a key role in their temporal evolution. 


IV. PHASE VELOCITY TRANSIENT 
DYNAMICS 

Many studies on shear flows have shown the impor¬ 
tance of the early time dynamics, that in principle can 
lead to large transient energy growth long before the ex¬ 
ponential mode becomes dominant [13, 41-46]. Transient 
dynamics offers a variety of different behaviors and phe¬ 
nomena which are not easy to predict a priori. In the lit¬ 
erature transient dynamics is generally analyzed in terms 
of energy amplification factor, while lesser attention has 
been paid to the phase speed evolution. However, the 
temporal evolution of the phase speed can lead to inter¬ 
esting considerations and useful information. In partic¬ 
ular, we can understand why and when different time 
scales appear inside the transient (see figures 4-7 below) 
and how these scales are linked to the features of the 
eigenvalue spectrum (for instance the spectral width). 
In fact, in the following, it can be observed that the 
sequence of different transient phases leading to the fi¬ 
nal temporal asymptote can be better observed through 
the phase speed rather than through the kinetic energy 
growth factor. 

In the context of our formulation, proper initial condi¬ 
tions on V and Coy have to be associated to the system of 
equations 11-12. We have computed the phase velocity 
through the transversal velocity component v, see Eqs. 
19-21 but we could have considered any component of 
velocity and vorticity. Moreover Eq. 11 is homogeneous 
and it can be demonstrated that the eventual introduc¬ 
tion of an initial transversal vorticity does not actually 
affect the perturbation temporal evolution, see Ref. [25] . 
For this reason, the initial vorticity a)j,(0, y) is set to zero 
for all the simulations here considered. As for the velocity 
field, we consider four different initial conditions, v{0,y) 
that we selected on the basis of the symmetry and the 
location of the perturbation with respect to the higher 
or lower vorticity region in the basic flow. 

For a fixed wavenumber and Reynolds number, the 
combination of these two attributes produces a different 
level of excitation of the least stable subset of eigenvalues. 
In particular, for the base flows here considered, the ex¬ 
cited least stable eigenvalues can belong to both the left 
and the right branches of the spectrum, see panels (d,h) 


















FIG. 3. Initial velocity disturbances Do and base flows (dashed 
lines). The amplitudes have been scaled for clarity. Top pan¬ 
els: channel flow. Bottom panels: wake flow. The blue lines 
represent symmetric disturbances (right panels), while the an¬ 
tisymmetric are in yellow (left panels). Thick lines indicate 
central initial disturbances, thin curves represent the lateral 
ones. 

in Fig. 4. After a very first stage where the contribu¬ 
tion from the most damped eigenmodes disappears very 
quickly, and since the growth/decay rates in the subset 
of the least stable eigenvalues can be vey close, a second 
stage may set in where a sort of intermediate temporal 
asymptotics lasting hundreds or more time scales inter¬ 
venes. This intermediate term leads to the final state, 
where the phase velocity takes the value corresponding 
to the last symmetric (or antisymmetric) eigenvalue, de¬ 
pending on the parity of the initial condition, see Fig. 
4(d,h). 

The exact expressions of the initial conditions we use 
are given in Table III and they are shown for both base 
flows in Fig. 3. In synthesis, we considered four dif¬ 
ferent initial conditions for both the plane wake and 
channel flow. For the wake flow they will be named SI, 
AI, SO, AO (S=symmetric, A=antisymmetric, I=inside- 
wake, 0=out-of-wake, see bottom panels of Fig. 3). For 
the channel, which is bounded, we name them as SC, AC 
, SW, AW (C=central, W=wall, see top panels of Fig. 
3). We recall that the for the channel flow the high shear 
region is located near the wall while for the wake flow it 
is located in the central region of the open domain that 
computationally is typically from 20 to 100 times larger 
than the base flow. 

An overview of possible transients is presented in Fig. 
4 where the evolution of three-dimensional perturbations 
is in terms of amplification factor G, phase velocity c, and 
temporal growth rate r. The transient of perturbations 


TABLE III. Initial conditions imposed on the velocity. Do = 
v{0,y). For the channel flow they are named SC/AC: 
symmetric/ antisymmetric and central. SW/AW: symmet¬ 
ric/antisymmetric and wall-located (the initial condition has 
large variations close to the channel walls). For the wake flow, 
SI/AI: symmetric/antisymmetric and inside-wake. SO/AO: 
symmetric/antisymmetric and out-of-wake (see Fig. 3). 


I.C. 


Wake flow 

Channel flow 

SI/SC 


Do = cos{y) 

Do = e 0 ^ cos(3i/) 

AI/AC 


Do = sin(2/) 

Do = e~o^ sin(3y) 

SO/SW 

Do 

= g-G-io)^ _l_ ^-(y+iof 

Vo = (1 - 

AO/AW 

Do 

= g-(y-iof _ ^-(y+wf 

Vo = 1/(1 - 


with initial conditions as in Table III are shown for both 
the base flows and for k < kd- The results obtained by the 
IVP are compared with the eigenvalue spectra given by 
the modal analysis. The trends for small wave numbers 
are shown here, as they could be more easily observable in 
the laboratory. However, in the summary scheme of Fig. 
9 the transient behavior of disturbances with wavenum¬ 
bers larger than kd is shown. For the two kind of shear 
flows, the amplification factors and the temporal growth 
rates are reported in panels (b,c,f,g) of Fig. 4 . For the 
wake flow, we can observe that with an initial condition 
of the kind SI - symmetric and localized inside the base 
flow shear region - and AO - antisymmetric external to 
the base flow shear region - the long term exponential 
behavior is reached after few temporal scales. On the 
contrary, with the other kind of initial conditions - SO 
and AI - the transients can last up to many hundreds 
of time units. In these cases quite far along within the 
transient the phase velocity abruptly changes its value, 
see panels (e). For the wake, initial perturbations SO 
slow down, that is the phase velocity decreases after the 
jump. Perturbations AI, instead, experience an acceler¬ 
ation. 

The same behavior is observed in the case of the channel 
flow, but here the role of the initial conditions is reversed. 
In fact, SW and AC show an exponential trend after few 
time scales, while SC and AW have a long transient char¬ 
acterized by the presence of phase velocity jumps (see 
Fig. 4a). Also in this case, the phase velocity of the sym¬ 
metric perturbation jumps to a smaller value while in the 
antisymmetric case it shifts to an higher value. The dif¬ 
ferent behavior of the two base flows can be traced back 
to the structure of the eigenvalue spectra and the associ¬ 
ated eigenmodes. Indeed, it is known that for the channel 
flow the eigensolutions of the left branch (A-branch) are 
wall-modes while those of the right branch (P-branch) are 
central modes. The opposite is true for the wake flow, 
whose spectra are made of a discrete set (left branch) 
of inside-wake modes, and a continuous branch of out- 
of-wake modes. As for the channel case, both branches 
may contain symmetrical and antisymmetric eigenfunc¬ 
tions, see panels (d,h) of Fig. 4. 
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FIG. 4. Temporal evolution of perturbations in terms of the phase velocity c (a,e), the amplification factor G (b,f), and the 
temporal growth rate r (panels c,g), for different initial conditions. The eigenvalues spectra are shown in panels (d, h). Left 
column: channel flow with Re = 6000, ij) ~ '^/4 and fc = 1. Right column: wake flow with Re = 100, <j) = tt/G and k = 0.7. 
The phase velocity is computed at yo = 0.5 for the channel and at i/o = 1 for the wake flow (NB: the reference lengths are 
the channel half-width and the cylinder diameter). The quantity Tc, see (a) and (e), indicates the temporal periodicity of the 
phase velocity fluctuations observed in the early and intermediate term and corresponds to the ratio Tc = 27r/[CTr„,j„, — 
as shown in panels (d) and (h). 
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As already reported in [25], besides the two temporal 
scales associated to the value of the phase velocity before 
and after the jumps, we observe a further periodicity, T^, 
related to the temporal modulation of the phase velocity 
during the early and intermediate terms, see Fig. 4(b,f). 
In this study we tried to relate this periodicity with the 
structure of the eigenvalue spectrum and verified that 
it is inversely proportional to the width of the range of 
frequencies given by the modal theory, Tc = 27r/[cri.^^^ — 

The system presents other two temporal scales: the 
external scale related to the base flow and the length of 
the transient (which can be determined by observing the 
time instant beyond which the growth rate r, and the 
angular frequency w, are both constant). Therefore, for 
a given wavenumber, from the phase speed transients it 
is possible to observe five different time scales. 

Fig. 4 shows that symmetric and antisymmetric pertur¬ 
bations do not have the same temporal asymptotic be¬ 
havior. More precisely, symmetric perturbations are al¬ 
ways less or equally stable than the antisymmetric ones. 
This means that if one considers mixed initial condi¬ 
tions, the symmetrical component will always prevail. Of 
course, in laboratory experiments, one cannot assume the 
precise symmetry because residual disturbance cannot be 
completely suppressed. As a consequence, the phase ve¬ 
locity that can be eventually observed is that given by the 
symmetric part of the initial condition. In an attempt to 
reproduce what might happen in the laboratory, one can 
build up a mixed initial condition with prevailing anti¬ 
symmetric component, by adding to the AI and the AO 
conditions a small noise which models the presence of a 
small symmetric part. The evolution of such a distur¬ 
bance is shown in Fig. 5: after an early transient, the 
phase velocity experiences a jump (as it happens for AI 
and AO) and then the solution keeps the antisymmetric 
shape until the symmetric part becomes dominant. The 
temporal asymptote is eventually reached (t < 2900) and 
it is announced by a second frequency jump. In this case 
the solution is a symmetric growing mode. 


V. NEAR-SIMILARITY OF VELOCITY 
PERTURBATIONS 

The observations of frequency jumps yield an interest¬ 
ing result: the perturbations temporal evolution have a 
three-part structure, with an early stage, an intermediate 
stage and an asymptotic stage. This multistage structure 
and the nature of the phase velocity jumps discussed in 
the previous section can be better highlighted by consid¬ 
ering the temporal evolution of the velocity modulus jDl 
normalized to its instantaneous oo-norm, see Figs. 6 and 
7. 

The early transient is the stage where the perturba¬ 
tion is most affected by the fine details of the initial 
condition and represents a period of adjustment. Dur¬ 
ing the intermediate stage the perturbation evolves al¬ 



I ’ 0 __,_iZ_ 

1000 ^ 2000 3000 



FIG. 5. Temporal evolution of mixed symmetric- 
antysimmetric initial conditions for the wake flow with Re = 
100, k = 0.4, 0 = 7r/3 and xq = 50. A small uniform random 
noise of magnitude 0(10“^) is added to the initial conditions 
AI and AO in order to simulate the transient life of a gen¬ 
eral initial disturbance, (a) Growth rate, (b) Phase velocity. 
The dashed line represents the asymptotic trend given by the 
least stable Orr-Sommerfeld eigenvalue. More information 
about the solution is shown in Fig. S3 in the SM online [37]. 
(c) Temporal evolution of the velocity profile modulus for the 
noisy AO initial condition. The profiles are normalized with 
respect to their peak value. 


most exponentially: the phase velocity takes the final 
constant value, the transverse velocity profiles maintain 
a near self-similar nature in time and the growth fac¬ 
tor changes very slowly. This stage appears as a kind of 
temporal intermediate asymptotics and can be in general 
considered extinguished only when both the frequency 
and the temporal growth rate become constant so that 
the long term regime is reached. This transient structure 
is common to all perturbations. In the case in which 
the frequency transient presents one or more jumps, the 
beginning of the intermediate asymptotics can be envis¬ 
aged where the near-selfsimilar features first appear. In 
particular, this first appearance can be observed when 
the initial conditions leads to two phase velocity jumps 
and in between these jumps the phase speed keeps a con- 















11 


slant value for many (possibly, up to hundreds or thou¬ 
sands) time scales. See Figs. 6 and 7 and the qualita¬ 
tive schematic shown in Fig. 9 where the different kinds 
of phase speed transients that we observed in this work 
are resumed. It should be noted that the establishment 
of self-similarity requires a constant wave propagation 
speed. 

In Figures 6 and 7, we present the cases of Re = 100, 
Xq = 50, k = 0.7, (j) = 7r/6 with SO and AI initial condi¬ 
tions for the wake flow and of Re = 6000, fc = 1, ^ = 7r/4 
with SC and AW initial condition for the channel flow. In 
the wake case the similarity is showed by normalizing the 
solution by the peak value of the profiles (the oo - norm) 
and by normalizing the lateral coordinate by the instan¬ 
taneous width of the perturbation profile (the distance 
from the axis where Ihl/UDHoo = 0.01). In the channel 
case, since the lateral diffusion is blocked by the walls, to 
highlight the similarity, it was sufficient to normalize the 
solution by the profile peak value. 

In the reported cases, the wake perturbation width 
scales in time as with p « 0.42. The general trend 
for p is shown in Fig. 8. One can observe that the expo¬ 
nent decreases as the Reynolds number and increases as 
the wave obliquity and the polar wavenumber. For very 
low values of the Reynolds number, when the inertial ef¬ 
fects become very little, the exponent is expected to take 
the value 0.5. Note that Re « 20 is the smallest reli¬ 
able value for which the wake flow can be represented in 
terms of a matched asymptotic expansions solution valid 
in both the spatial intermediate and far field, see Ref. 
[21] for major details on the wake base flow computa¬ 
tion. The dependence on the obliquity angle is weak at 
low angles, while the diffusive scaling occurs for per¬ 
turbations orthogonal to the base flow {(j) = 7r/2) since 
in this case the convective transport does not matter, 
see Eqs. 11-12. Moreover, as expected, we observe that 
short waves present a more diffusive behavior than long 
ones, see the panel (b) in Fig. 8 . Fig. 9 shows a qual¬ 
itative scheme of the phase velocity transient type, for 
different combinations of wavenumber - type of initial 
condition. The diagram distinguishes the cases where 
the perturbation is accelerated or decelerated after the 
jump. For instance, for the wake, in the case of expand¬ 
ing perturbation along the y direction, as shown in Fig. 
7 b, the change in the perturbation shape occurs when 
the phase speed starts to shift to the value specific of 
the least stable eigenvalue. Other examples of acceler¬ 
ation/retardation of the wave propagation can be find 
in Figs. 4, 5, 6, 7. It is also noteworthy that the time 
interval where the frequency/phase velocity jumps are 
observed is rather distant from the initial moment, and 
for this reason it could be observed in the laboratory. For 
instance, in the case of the wake flow with Re = 100 and 
AI type of initial condition (see panel e in Fig. 4), if 
we consider a cylinder with D = 2 cm, the jump arrives 
after nearly 4 minutes if the fluid is air, and after nearly 
50 minutes if the fluid considered is water. In the case 
of the channel flow, Re = 6000 and SC type of initial 



Channel flow, 

Re=6000, k= 1, (|)=jt/4 

y„=0.5 

AW 

- SC 
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FIG. 6. Channel flow, Re = 6000, fc = 1, 0 = Tr/fl, ini¬ 
tial conditions AW (yellow, thin curves) and SC (blue, thick 
curves), (a) Temporal evolution of the phase velocity, com¬ 
puted at yo = 0.5. A jump is observed at t ~ 100 for SC 
and at t ~ 600 for AW. (b) Evolution of the velocity profile 
modulus for the SC initial condition. The profiles are nor¬ 
malized with respect to their peak value, (c) Evolution of the 
normalized velocity profiles for the case of AW initial condi¬ 
tion. The phase velocity shifts indicate a substantial change 
in the velocity profiles. This transition can be more or less 
abrupt depending on the simulation parameters, and the re¬ 
gions before and after the jump can show a nearly self-similar 
behavior. 

condition (see panel a in Fig. 4), if we consider a chan¬ 
nel with h = 15 cm, the jump arrives after nearly 25 sec 
if the fluid is air, and after nearly 6 minutes if the fluid 
considered is water. 


VI. CONCLUDING REMARKS 

We observed the existence of a wavenumber threshold 
kd that splits the range where the least-damped waves 
in the long term disperse from the range where they are 
non-dispersive. The transition between dispersive and 
non-dispersive behavior is related to the fact that for 
k < kd the least-stable Orr-Sommerfeld eigenvalue be¬ 
longs to the left branch of the spectrum, while for k > kd, 
it belongs to the right branch. We considered the tem¬ 
poral kinetic energy evolution in the wavenumber space 
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FIG. 7. Wake flow, Re = 100, k = 0.7, (p = ^r/S, xq = 
50, initial conditions AI (yellow, thick curves) and SO (blue, 
thin curves), (a) Temporal evolution of the phase velocity, 
computed at yo = 1- A jump is observed at t « 400 for SO and 
at t « 800 for AI. (b) Evolution of the velocity profile modulus 
for SO. The profiles are normalized with respect to their peak 
value, (c) Evolution of the normalized velocity profiles for the 
case of AI ic. It can be noticed that the solution for the wake 
flow admits the presence of near-selfsimilar terms (see also 
Fig. 9) of expanding velocity profiles. It is also interesting to 
note that the rapid transition happens when the perturbation 
reaches the wake from the outer region (as the SO case, panel 
(b)), or viceversa as for the AI case (c). The perturbation 
width follows a power law p ~ 0.42 in the above cases. 


and distinguished the rate of change of kinetic energy due 
both to viscous and convection effects. It is observed that 
in the wavenumber space the kinetic energy shows an ex¬ 
change of mode aX k = k^. In fact, the kinetic energy of 
the solution follows that of the least stable mode which 
is located on the tip of the left branch of the spectrum 
for k < kd and is located on the tip of the right branch 
for k > kd- The transition in the dispersion relation also 
helps to explain some features of the transient dynam¬ 
ics. Indeed, initial conditions with a wave number higher 
or lower than kd have different transients. The possible 
types of phase velocity transients that we observed in this 
study can be classified into four categories: 

1. short transients where the phase velocity reaches 
its long-term asymptote smoothly (no jumps) 

2. long transients where the phase velocity jumps to 



FIG. 8. Scaling exponent p for the temporal evolution of the 
wake perturbation width in the intermediate term at xq = 50 
and with AI initial condition, (a) The exponent is represented 
as a function of the Reynolds number and it is parameterized 
with the obliquity angle in the range [0,7r/2]. (b) The depen¬ 
dence on the obliquity angle is shown for Re = {30, 70,100} 
and k = {0.7,1.1}. 


a higher value and the wave propagation is acceler¬ 
ated 

3. long transient where the phase velocity jumps to a 
lower value and the wave propagation is decelerated 

4. long transients where more than one jump occurs 
during the phase velocity transient. 

Jumps are in general accompanied by a modulation of 
the phase velocity which is observed both before and af¬ 
ter the instant when they take place. This modulation 
shows a temporal periodicity that is highly correlated 
with the width of the range of frequencies given by the 
modal theory. The type of transient resulting of course 
depends also on the imposed initial condition. Initial 
conditions can be classified into four types according to 
their symmetry and to the centrality or laterality of the 
regions where the initial condition takes its largest vari¬ 
ation. The observation of phase velocity transients leads 
to identifying three distinct stages in the time evolution 
of the perturbations. There is always a first part, the 
early transient, which is heavily dependent on the initial 
condition, a second much longer part, the intermediate 
transient and a third part, the long-term state that is 
reached when both the phase velocity and the growth 
rate take their final constant values. Inside the intermedi¬ 
ate stage, interesting and non-intuitive, the perturbation 
evolves almost exponentially: the phase velocity takes a 
constant value, the transverse velocity profiles maintain 
a near self-similar nature over time and the growth factor 
changes very slowly. In the wake flow, which is a system 
slowly evolving in space, our simulations highlight a lat¬ 
eral growth of the perturbation profiles which follows a 
temporal power law scaling where the exponent changes 
not only with the Reynolds number, but also with the 
wavelength and wave-angle. In the case of purely sym¬ 
metric or antisymmetric initial condition, the beginning 
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P base velocity 
qualitative transient 


Wake flow 


Channel flow 


No jumps 
Short transient 


Initial conditions SI or AO, 
andk<k, 

(non-dispersiveness)l III 1111 
Initial conditions SO or AO, 
and k > k^ 

(dispersiveness) mrff 


Initial conditions SW or AC, 
1 


and k < k, 

d 

(non-dispersiveness) ^ 

Initial conditions SC or AC, 
and k > k, HL iiii 
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dispersiveness 


Initial conditions Al, 
and k < k^ * 

Initial conditions SI or AI, 
and k > k, 


Initial conditions AW, 
and k < k, * 


Initial conditions SW or AW, 
and k > k, 
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Initial conditions SO, 
and k < k. 
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Initial conditions SC, 
and k < kj 
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decelerating, dispersiveness 



Mixed initial conditions 
with dominant antisymmetric 
component, and k < k^ 


ffF 


Mixed initial conditions 
with dominant antisymmetric 
component, and k < k^ 


FIG. 9. Qualitative scheme for the phase velocity transients. Left column-, type of transient; the black lines represent the 
phase velocity evolution over time (refer to Figs 4-7 in the paper), while the red sketches indicate the presence of self-similar 
perturbation profiles, in the intermediate term, for the wake case (please refer to Fig 7). Mid column-, wake flow; in each 
cell the conditions under which the corresponding kind of transient realizes are indicated, in terms of initial condition type 
(SI=symmetric in-wake, AI=antisymmetric in wake, SO=symmetric out-of-wake, AO=antisymmetric out-of-wake) and the 
wavenumber value, namely k < kd or k > kd- Right column-, channel flow. (SC=symmetric central, AC=antisymmetric central, 
SW=symmetric lateral, AW = antisymmetric lateral). With the blue sketches we represent, for each case, the shape of the 
perturbation in the long term (temporal asymptotic conditions). The green star is to remind that, even for k < kd, the least 
stable antisymmetric mode belongs to the right branch of the spectrum and it is non-dispersive. 


of the intermediate transient is marked by the phase ve¬ 
locity jumps, if there are any. Transients show two phase 
speed jumps when the initial condition is antisymmetric 
and noisy. In this case, near self similarity begins beyond 
the first jump where the propagation accelerates. 
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Appendix A: A Chandrasekhar eigenfunction 
expansion method to solve the linearized initial 
value problem 

In this section, based on an eigenfunctions expansion 
in terms of the Chandrasekhar functions [47] , a Galerkin 
method to solve the Orr-Sommerfeld and Squire initial 
value problem is presented. The method is valid for ev¬ 
ery parallel base flow with homogeneous boundary con¬ 
ditions and is conceptually an extension, specihc to non 
modal calculations, of the use of Chandrasekhar’s func¬ 
tions, which satisfy both the vanishing of the normal ve¬ 
locity and its normal derivative on wall boundary con- 
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ditions of viscous flows. In particular, it is a three di¬ 
mensional extension of the method by Gallagher & Mer¬ 
cer [48] which was based on the use of Chandrasekhar’s 
functions to solve the bi-dimensional modal perturbative 
problem for the Couette flow. 

As a starting point, the IVP in the wall-normal velocity 
and vorticity form is considered (this is equivalent to the 
system 11-12). 


dtdyV — k^dtv -\- iaU{y)dyV — iak^U{y)v — iaU"{y)v 

(Al) 


dt(2}y + iaU {y)u}y 


1 

Re 



—ijU (yYv 
(A2) 


^{y = ±1: i) = dyviy = ±1, t) = Uy{y = ±1, t) = 0 

(A3) 

v{y,t = 0) =vo{y) Wy(y, t = 0) = a)y„(y) (A4) 


1. Solution to V equation 


The solution of (Al) can be expressed as a generalized 
Fourier expansion with time-dependent coefficients: 


'Ky,t) = '^Cn{t)Xn{y) ye [-1,1], 


(A5) 


n=l 


where Xn{y) are orthogonal functions, and the following 
inverse transform applies: 


!^iv{y,t)Xn{y)dy 

Cn(c) = -• 

/_iX„(y)X„(y)dy 


(A6) 


Since in the initial value problem both the initial condi¬ 
tion and the boundary conditions need to be imposed, 
it is worthwhile to consider functions that satisfy the 
boundary conditions. Moreover note that the coefficients 
Cn of the series are in general complex, since v is complex¬ 
valued and the spatial modes are considered as real. The 
particular orthogonal functions which we use are those 
defined by the following fourth order eigenvalue problem 

^^X{y) = \^X{y) ye [-1,1] (A7) 

X(y = ±l)=0 ;^X(y = ±l)=0. (A8) 

dy 

Two different sets of eigenvalues and the correspond¬ 
ing eigenfunctions are found, respectively odd and even, 
by numerically solving the following transcendental equa¬ 
tions 


tan{\n) — tanh{\n) = 0 {odd set) (A9) 
tan{\n) +tanh{Xn) = 0 {even set). (AlO) 




FIG. 10. The basis eigenfunctions 


The corresponding normalized eigenfunctions (figure 10) 
are 


sinh{X„y) 
sinh{Xn) 

{odd set) 




sin{Xny) 

sin{Xn) 


n = l,3,5..,fV- 1 
(All) 


COsh{Xny) 
COSh{Xn) 

{even set). 



cos{Xny) 

COs{Xn) 


n = 2,4,6..,Af 

(A12) 


Similar functions, in a different domain, have been used 
in the study of the circular Couette flow between coaxial 
cylinders [47, appendix V]. 

Since the imaginary and the real part of the solution v 
usually have opposite parity, independently on the initial 
condition, both the odd and the even sets are necessary 
to completely describe the problem and obtain the cor¬ 
rect result. 

In the following paragraphs a compact notation for the 
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space derivatives is introduced. In order to simplify the 
reading, the y-derivatives will be indicated with a sub¬ 
script. The temporal derivatives will be indicated explic¬ 
itly or with a dot. 

The numerical solution to the i) equation (Al) is 
obtained by applying the variational Galerkin method. 
Truncating the series (AS) at N functions and substitut¬ 
ing yields 


Dr. 


= {X„,Xm)=Sr 


(A17) 


— {^riyy ) Xm) — 


(A18) 


{ I A M 

[An^n 

0 


^ml^n 


if (n -|- m) is even, n ^ m 
if (n + in) is odd 
if n = TO 


a, 


JV d ^ 

7) = H ^Cn{t)X„ 


dt ' ' ^ dt 

n—1 n=l 

N 


N 


+ iaU{y) YyCu{t)Xr.yy-^akMy)Y.^r.{t)X. 

L 

N I ^ 

Cn{t)Xn — — 5]] Cn{t)X.„yyyy 
n—1 


. d^U{y) 

— la 


dy2 


n—1 


.,n ^n^TTi^n 

(AI9) 


(A20) 

= {U{y)X^,Xm) 

(A21) 

[/(3) _ /*^ X X ) 

(A22) 



N 

C'n(t)Xnyy 

n—1 


Re 


N 

^ ^ C'n(j^)^n 

n—1 


(A13) 


where 

cosh{2X„) — cos(2An) 
sinh{2\n) — sin{2Xn) 


lim /j,„ = 1. (A23) 


The error functional e is minimized when it is orthogo¬ 
nal to the space of the linearly independent trial functions 
X„ with n = 1,2, ..N. In this context, given two func¬ 
tions u{y) and v{y) with y G 11 = [—1, 1], the following 
definition of scalar product applies 

{u,v) = / u-vdy. (A14) 

Jq 


It is convenient to express the ODE system (AI6) in 
a compact notation: in the following, vectors will be in¬ 
dicated either explicitly using braces or with bold lower 
case letters; matrices will be indicated with bold capital 
letters; constants with roman capital letters and physical 
parameters in italic. The system can be written as 

Hc-Gc = 0, (A24) 


With above notation the Galerkin orthogonality condi¬ 
tion is expressed as 

{e,Xm) = 0 to = 1,2,...,W (A15) 

The Orr-Sommerfeld PDE is now reduced to a system 
of N ODEs of the first order, where the time dependent 
coefficients c„(t) are the only unknowns. 


N , N 

0 = y^ ~^C,n{t){Xnyy, Xra) — k 5^ 

n—1 n—1 


_d 

dt 


(t) {Xji , Xm) 


where 


H = S - D 


G = —ia U(D -I- iak'^ XjG) 

+ ia U(3) 


-t- 



2fc2 

Re 


Re 


D, 


(A25) 

(A26) 


where D = [Dm,n] etc., i.e. the element Dm,n is placed 
at the column and at the to*^ row of the matrix. H 
is invertible, so denoting A = G yields 


c-Ac = 0. (A27) 


N 


N 


-b ia 


n—1 

N 


n—1 

N 




n—1 
l2 N 


n—1 


N 


2k‘^ 

A 'y ^ C„(t) {X„yy , Xm) — y ^ C„(t)(A„, Xm) 


Re 


The complex eigenvalues Ci of A constitute the spec- 
);rum of the Orr-Sommerfeld equation, and the analytic 
solution to eq. A27 can be easily implemented and com¬ 
puted numerically(the QR method is used to compute 
the eigenvalues, in the matlab environment). 


2. Solution to the nonhomogeneous ujy equation 


n, TO = 1,2,3,..., A. (A16) 

The scalar products can be evaluated analytically or 
computed by numerical integration: 


In order to solve the Squire IVP, a set of normal func¬ 
tions different from the one adopted for the velocity is 
needed, since the second order PDE only requires uiy to 
vanish at the boundaries, but not its first derivative. A 
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simple choice for the basis functions, here adopted, is the 
following 


Yn = sin{^ny) 
Yn = COsi^nV) 


where 

u = 

^n = 


{n + l)7r 
2 

(n — l)7r 


n = 1, 3, 5, ...N — 1 {odd set) 

(A28) 

n = 2,4, 6, ...A^ {even set), 

(A29) 


n = 1, 3,5, ...N — 1 {odd set) 

(A30) 


n = 2,4,6,...A^ 


{even set). 
(A31) 


Also in this case, note that two sets of eigenfunctions are 
put together to form a unique set, since both are neces¬ 
sary to completely describe the complex-valued normal 
vorticity. The general solution is then obtained as the 
sum of a particular solution and the solution to the 
corresponding homogeneous equation ujyf^ 

<^v{y,t) = <^yh.{y,t) + <^yp{y,t) (A32) 


<^yiy,t) = + bpJ{t)Y„{y). (A33) 

n—1 

Applying the Galerkin method yields to the following 
forced ODE set 


b-( -jaU* + —S*-^D* 
Re Re 


b-G*b = Bc, 



(A34) 

(A35) 


where the eigenvalues of G* constitute the spectrum 
of the Squire equation and 


D*^^^ = {Y^,Y^)=5m,n (A36) 

S*^,n = {Ynpp,Y^)=-ej^,n (A37) 

U*^^^ = {U{y)Y^,Yrr.) (A38) 

Fkn = (^^Xr„Y^). (A39) 


The homogeneous solution bh of (A35) is analytically 
known, while a particular solution bp of the following 


form is sought 

N 

bp„{t) = X! (A40) 

i=i 

where Unj are constants and Oj are the eigenvalues of G, 
in such a way the forced solution has the same spectral 
content of forcing term. The coefficients are obtained 
through the solution of N algebraic sets, after the com¬ 
putation of c(0). 

The Galerkin method was first applied to the Orr- 
Sommerfeld modal equation by [49]. They used nor¬ 
mal functions that guarantee a 1/A^^ convergence ra¬ 
tio. Gallagher [48] used, for the modal problem, the 
Chandrasekhar-Reid functions and the error decreased 
as 1/A^^ with N ^ oo as shown in [50]. The fifth order 
of accuracy is ensured for the present formulation as well, 
as shown in figure 11. 

The method results to be fast and accurate in time 
and space. Since the time evolution is analytically rep¬ 
resented, the complete wave transient, up to the asymp¬ 
tote, can be simulated without typical drawbacks of time¬ 
marching techniques. Arbitrary initial conditions can be 
specified for bounded flows. The limits of the method are 
related to the non-normality of the Orr-Sommerfeld and 
Squire operators. The non-normal effects act on the nu¬ 
merical procedure by worsening the condition number of 
the eigenvector matrices. Anyway, the sensibility of the 
spectrum (especially at high Re and k values) is a prop¬ 
erty of the stability operator and it is thus independent 
on the numerical scheme. 


Appendix B: Spectra computation 

Two different numerical methods are used in the 
present work to compute the spectra of the Orr- 
Sommerfeld equation: - the fourth order hnite differences 
collocation scheme [51] and - the present non modal three 
dimensional version of the Gallagher and Mercer (1962) 
method which is described here in Appendix A. Here, we 
report a comparison with literature results as a validation 
of our spectral calculations, see Fig.12. For unbounded 
flows it has been shown by Grosch [2] that a continu¬ 
ous spectrum can be analytically found, if the bound¬ 
ary conditions are relaxed to v bounded as ?/ —>■ oo. If 
finite-norm boundary conditions are imposed and there¬ 
fore only the class of decaying solutions as \y\ —>■ oo is 
considered, the continuous part of the spectrum is ap¬ 
proximated in a discrete way. If the boundary condi¬ 
tions are imposed far from the wake, the approximation 
is very good. The Galerkin method with Ghandrasekhar 
functions described above was successfully adapted to the 
wake flow and to the boundary layer flow. Since no spec¬ 
tra with our wake base flow (see [21]) have been found in 
literature, the schemes have been validated with the Bla- 
sius boundary layer flow (see figure 12, panel a). Eventu¬ 
ally, for the channel flow, the comparison with a hybrid 
spectral collocation method based on Ghebyshev polyno¬ 
mials [8, A.6] is shown in figure 12, panel b. 
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FIG. 11. Maximum and rms of the absolute error of v 
(left panel) and Wy (right panel) as a function of the num¬ 
ber of modes N for channel flow with to = 100, Re = 
1000, k = 2, (j) = 80° and symmetrical initial condi¬ 
tion. Continuous line: real part. Dashed line: imagi¬ 
nary part. Magenta line: maximum absolute error. Blue 
line: rms of the absolute error. Black line: accuracy 
trend N~^ [50]. Since the exact solution is not known, 
the residuals are defined as the difference between the so¬ 
lution and an accurate solution computed with 350 modes. 

£a{y, t) = \vN{y ,t) - VN=35o{y,t)\, rms{ea){t) = 

£l{y,t), max{ea){t) =raaxy.{ea{y,t)) 
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FIG. 12. Eigenvalues spectra {a = ar + ioi) of the Orr- 
Sommerfeld equation. Comparison of different numerical 
methods: a 4*^ order finite differences scheme on uniform 
grid (blue circles); a 5^^ order Galerkin method on nonuni¬ 
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